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A stochastic nonlinear partial differential equation is built for two different models 
exhibiting self-organized criticality, the Bak, Tang, and Wiesenfeld (BTW) sandpile 
model and the Zhang's model. The dynamic renormalization group (DRG) enables to 
compute the critical exponents. However, the nontrivial stable fixed point of the DRG 

■ transformation is unreachable for the original parameters of the models. We introduce 
an alternative regularization of the step function involved in the threshold condition, 
which breaks the symmetry of the BTW model. Although the symmetry properties of 

■ the two models are different, it is shown that they both belong to the same universality 
class. In this case the DRG procedure leads to a symmetric behavior for both models, 
restoring the broken symmetry, and makes accessible the nontrivial fixed point. This 
technique could also be applied to other problems with threshold dynamics. 
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In the last decade much attention has been paid to the phenomenon known as self-organized critical- 
ity (SOC). Bak, Tang, and Wiesenfeld (BTW) jjj studied a cellular automaton model as a paradigm 
for the explanation of two widely occurring phenomena in nature: 1// noise and fractal structures. 
Both have in common a lack of characteristic scales. The SOC models, although not always show 
1// noise, have no characteristic scales too, and this scale invariance suggests that these systems are 
critical in analogy with classical equilibrium critical phenomena, but in SOC one deals with dynami- 
cal nonequilibrium statistical properties. Moreover, the system evolves naturally to the critical state 
Sh , without any tuning of external parameters, that is, in a self-organized process. 

Several cellular automata and coupled map lattices models exhibiting SOC have been reported in 
the literature. In the original sandpile model of Bak et al. [0 the system is perturbed externally by 
a random addition of sand grains. Once the slope between neighboring cells has reached a threshold 
value, sand is transferred between them in a fixed amount. Taking this model as a reference, different 
' dynamical rules have been investigated leading to a wide variety of universality classes. Continuous 

variables with a full transfer from a cell instead of a fixed discrete amount ||-[|, directed flows Q, 
threshold condition imposed on the height, on the gradient, or even on the laplacian 0, anisotropy 
are a few examples. These randomly driven models do not exhibit SOC when the interaction 
rules are not conservative H. Later on, other deterministically driven models have been introduced 
where conservation is not a necessary condition |10-15|. Much more recently, sandpile models with 
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deterministic perturbations but intrinsic randomness in the threshold dynamics have been used to 
reproduce experimental results on transport properties on rice piles |Tj| . The close connection between 
these sandpile models and interface depinning has been establish in Ref. ]l7| . 

Some authors have attempted to connect the randomly driven models to stochastic differential 
equations |l8| , |ll| . These continuous descriptions are built according to the symmetry rules obeyed by 
the discrete models in order to achieve a generic scale invariance condition p0| ]. Nevertheless none 
of them neither explicitly nor implicitly includes the threshold condition, which is one of the main 
characteristics of SOC models. On the other hand anomalous diffusion equations with singularities 
in the diffusion coefficient have been considered in order to study the deterministic dynamics of the 
avalanches generated in the critical state (2^j2^| . Latterly, a different approach has been introduced by 
Pictronero and co-workers, using a real-space renormalization procedure to determine the dynamical 
exponent as well as the avalanche size exponent. |B3|. 
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In a previous paper |2j| one of us studied two nonlinear stochastic differential equations derived from 
the discrete dynamical rules of two models with different symmetry properties. In principle one would 
expect, due to this reason, different critical behavior. However, it was shown analytically, by means of 
the dynamic renormalization group (DRG) |25|-|27j , that both models belong to the same universality 
class. The threshold condition was kept but the step function was regularized in order to allow a power 
series expansion. In the limit that recovers the threshold it was shown that the coupling constants 
that distinguish both models become decoupled of the common coupling constants; since the critical 
exponents only depend on the latter constants, one obtains the same values for both models. Once 
this equivalence was established, the most symmetric model was considered, showing that an infinite 
number of coupling constants was relevant below the upper critical dimension d c — 4; by expanding 
in the number of nonlinearities, the DRG procedure, up to first order in e = 4 — d, gave an estimate 
of the dynamical exponent close to the value obtained by scaling arguments and in the numerical 
simulations ^,|| . This value of the dynamical exponent is obtained when the flow in parameter space 
reaches the nontrivial stable fixed point; nevertheless, when taking into account the physical values of 
the parameters, they do not lie in the basin of attraction of the fixed point, thus making this computed 
value in some sense speculative, since it cannot be ensured that the flow in parameter space will be 
able to reach the attractor. 

Our goal in this paper is to complement the previous work in order to check the validity in the 
calculation of the critical exponents at the stable fixed points and to analyze the role played by 
symmetries in randomly driven SOC models and, in general, in other models where an infinite hierarchy 
of nonlinear terms is required. Our procedure also illustrates the effect of symmetry breaking in DRG 
calculations as a mechanism to make the attractors in parameter space accessible for the physical 
values of the parameters in the original equations. The continuum equation for the BTW and Zhang's 
models is introduced in Sect. II, as well as the alternative regularization which breaks the symmetry 
that distinguish both models. In Sect. Ill we develop the DRG procedure and show that this symmetry 
is irrelevant, in view of the fact that the nontrivial fixed point is not modified by this new approach. 
Moreover, the effect of symmetry breaking allows the flow of the original parameters to reach the 
nontrivial fixed point, where the critical exponents can now be computed. Finally we present our 
conclusions in Sect. IV. 

II. MODELS AND SYMMETRIES 

First, we describe briefly the dynamics of the two SOC models under consideration. The first model 
was originally proposed by Zhang Q), and it consists of a d-dimensional lattice in which any site 
can store some continuously distributed variable E. This variable, that we will call energy, can have 
different physical interpretations Q . The system is perturbed by adding a random amount of energy 
SE > at a randomly chosen site. Once a site reaches a value of the energy greater than the threshold 
value E c , this site becomes active, and transfers all its energy to its nearest neighbors. At this point 
the input of energy from the outside is turned off. The energy transferred to the neighboring sites 
can make them active, giving rise to an activation cluster or avalanche, that ends when all the sites 
have reached a value of the energy smaller than E c . It is only when the avalanche has stopped that 
energy is added again, otherwise the system remains quiescent. In this way there is a clear time-scale 
separation in the dynamics. The external noise acts in a slow time scale, whereas the avalanches evolve 
infinitely fast, in comparison. The second model differs from Zhang's model only in the amount of 
energy an active site transfers to its neighbors, which is a fixed amount E c , instead of its whole energy 
E. Therefore it is closer to the original sandpile model of Bak et al. but continuous in E. When 
5E is not random but fixed this difference becomes irrelevant. For this reason, it will be referred as 
BTW model. Notice that both models are conservative, in the sense that the added energy (always 
positive) is only dissipated at the open boundaries. 

The microscopic evolution rules can be written from time t (in the fast time scale) to t + 1 for each 
site i as 

Ei(t + 1) = Ei(t) - Ei(t)Q {Ei(t) -E c ) + -Y] E nn (t)G (E nn {t) - E c ) + &(t) (la) 

nn 

and 
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Ei(t + 1) = Ei(t) - E C Q (E t (t) -E c ) + -J2 E C Q (E nn (t) - E c ) + &(t), (lb) 

7171 

for Zhang's and BTW model, respectively. The sum runs over the q nearest neighbors of site i, 
labeled as nn, and the threshold condition enters through the Heaviside step function 0, defined as 
9(x < 0) = and <d(x > 0) = 1. Due to the continuous nature of the models the value <d(x — 0) is 
irrelevant and we can keep it undefined, by now. For the external noise that drives the system 

only when there are no active sites, one can formally write 

&(*) = SE5 iMt) J] [1 - E e )] , (2) 

where 8i tn u\ is the Kronecker delta symbol and n(t) is a random vector pointing the site of the lattice 
that will be perturbed with a random amount of energy SE (in the original BTW model 5E = E c /q). 
The product runs over all the lattice sites. 

When applying the DRG one deals with infinite systems and then the important effect of dissipation 
at the open boundaries is not taken into account. However in SOC models a distribution of absorbing 
defects through the lattice plays the same role as the open (absorbing) boundaries p8|, as we have 



verified through computers simulations 29 1 . Then, we can redefine our models in an infinite lattice but 
with a quenched distribution of defects. The results are not modified with this assumption. Another 
possibility is to consider that each site of an infinite lattice has a small probability of dissipating an 
amount of energy E c /q when it topples, instead of transferring it to a certain neighbor. This procedure, 
that represents the assumption of random boundaries, implies that when a site receives a toppling from 
some neighbor, it has a small probability of not to accept the amount of energy E c /q, that is lost fiof . 
This dissipation can be included as a new term in the noise, and Eq. (Q) has to be replaced by 

&(t) - SE5 iMt) H [1 6 (Ej(t) - E c )} - ]T ( m e(E m (t) - E c ), (3) 

Vj nn 

where £„„ is a dichotomous noise, taking value with large probability and value E c /q (dissipation) 
with a small one. If £„„ depends on t, i.e., ( nn = Cnn{i) we are dealing with annealed random 
boundaries, whereas if it only depends on the position, we have quenched random boundaries or 
absorbing defects. 

In terms of a rescaled energy E — E c — > E, and introducing a parameter Z to unify the description, 
we have for both models 

Ei{t+1)-Ei{t) = - {[ZE nn (t) + E c ] 9 (E nn (t)) - [ZEi{t) + E c ] 6 (E t (t))} + (4) 
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where Z = 1 for Zhang's model and Z = for BTW. Equation (|J) defines a stochastic coupled map 
lattice. Moreover, notice that the deterministic BTW equation displays invariance under a parity 
transformation of the order parameter, E — > —E. This one is the only symmetry that the BTW model 
does not share with Zhang's model. The common symmetries are invariance under spatial translations, 
rotations, and reflections, as well as conservation of the order parameter. 

Equation (Q) can be coarse-grained in order to obtain a continuum equation for the effective E(r, t). 
Then, by using the prescriptions for the temporal derivative and for the Laplace operator: 

8E(r,t) 



dt 



= aV 2 {[ZE(f, t) + E c ]Q{E{r, t))} + r,{r, t), (5) 



where a is a coefficient that depends on the lattice spacing, the unit time step and the coordination 
number q. The noise n(r, t) accounts for the effective external noise as well as for the internal noise 
that appears due to the elimination of microscopic degrees of freedom. 

Up to now, Eq. (^) describes truly the coarse-grained evolution of the system, but we have not yet 
characterized the noise n(f,t), which derives from ^(t). The product in Eq. (|^) makes n(r,t) to be 
a multiplicative noise that depends on the whole lattice state, and the problem is intractable. We 
are going to ignore the restrictions imposed by the step functions in Eq. (|^), breaking the time scale 
separation. Then the noise n(f, t) acts continuously in time and can provoke avalanches to overlap. 
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However for small enough noise this is very unlikely, and one can still identify avalanches in computer 
simulations. Moreover the dynamical exponent does not change with this assumption @], due to the 
fact that the added noise is orders of magnitude smaller than the energy transferred by the avalanche 
and then its dynamics is not affected. In this case we still have two time scales, although they overlap. 
In what follows ri(f,t) will be considered as an additive random process, including two effects, the 
external driving, always positive, and the dissipation at the (random) boundaries, always negative. In 
the statistical stationary state the random input of energy must equal on average the output at the 
boundaries. Then we assume that 

<T?(f,i) >=0. (6) 



In fact this is the same assumption done in all the studies of SOC by means of DRG [ |18|JlSH , and it is 
somehow equivalent to the stationary condition used in Ref. p3j . 

Moreover we are mainly interested in the spatio-temporal propagation of a perturbation through the 
system, that is, in measuring the value of the dynamical exponent. For this purpose we have to look 
at the system in the fast time scale, i.e., the scale of the evolution of the avalanches. In Ref. || it was 
argued that in this case one can understand the noise as a quenched Gaussian process uncorrelated in 
space, and then its correlation function is given by 

< 77(f,t)77(r',t') >= 2T5 d (r-r'). (7) 

When looking at the system at the slow time scale one cannot use this prescription for the noise, which 
has to be uncorrelated in time too, i.e., < r/(r, t)rj(r' , t') >— 2TS d (r — r')S(t — t'), and this prescription 
is mainly related to the interface roughness between avalanches |3l|] . 

Equations (||), (JsJ) , and (Q), together with the fact that the noise is a Gaussian process, define 
completely our model. However, the presence of the step function in Eq. (^) gives rise to a strong 
nonlinearity. A perturbative expansion of this equation can be performed if one regularizes the step 
function as 

0(E) = km f((3E), (8) 

and makes a series expansion of f{j3E) in powers of E ]22| , |24| ]. The function f(x) must be monotonously 
increasing with /(— oo) = and /(oo) = 1. Moreover we choose f(x) — 1/2 as an odd function, so 
f(0) = 1/2. Several functions of this type have been used in the literature, but that coming from the 
error function as 

1 + erf(x) 1 r a 
/ W = o = -~T I e v dy, (9) 



la 



allows a power expansion that has an infinite radius of convergence, in contrast with previous choices 
| f22||2~|l . In any case, the relevant results do not depend on the particular form of f(x). 

The regularization given by Eq. (|^) keeps the symmetry of the step function and therefore the 
invariance under a parity transformation in the BTW model. As an alternative regularization that 
breaks this invariance we propose the following: 

@(E) = lim f((3E + K) (10) 

(3— »oo 

with K an arbitrary constant. Although in the limit (3 — > oo we recover the step function, we do not 
recover its symmetry anymore, because <d(E = 0) = f(K) ^ 1/2 if K ^ 0, and this is the reason for the 
breaking of the symmetry in the BTW model. Now we perform a series expansion of the regularizing 
function f{(3E) in powers of E, obtaining 



9(E) = Urn J2a n (i3,K)E n , (11) 

71=0 

where the coefficients a n (P,K) are given by 



a n {(3,K) = j , (12) 
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being f^ n \K) the n-th order derivative of f(x) at x = K. Substituting the expansion ( |ll| ) into Eq. 
(||) we can write 

= DW 2 E{^t) + Y j \ n W 2 E n {f,t) +r,(f f ,t), (13) 

n—2 

where the effective diffusion constant D and the coupling constants X n (that make the equation non- 
linear) take different values depending on the model: 

D= lim a(E c fV(K)(3+Zf(K)), (14) 
A n = lim ^ (ejW(K) + Z Uf iK) ) , n = 2,3,...,oo. (15) 
On the one hand, for K — 0, since all the even derivatives ver 

ify /( 2ll+2 ) (0) = 0, all even coupling- 
constants vanish for the BTW model, whereas they do not for the Zhang's one. Using Eq. (jl^) this 
fact allows to verify the symmetry of the BTW model under the parity transformation of the order 
parameter E. On the other hand, for K ^ 0, the even coupling constants do not vanish in any case, 
and this fact constitutes the symmetry breaking for the BTW model. Then, under this condition, the 
only difference between both models is that the constants depend on f3 in a different way; however it 
is easy to see that in the limit j3 — ► oo both sets of constants are identical and then the Zhang's model 
and the broken-symmetry BTW model have to belong to the same universality class. This fact can 
only be shown for K =^ 0. Nevertheless, the fact of considering K = only introduces a difference in 
the value of O(0), which is irrelevant in a continuous model, and then one can include the (symmetric) 
BTW model in this universality class too. 

At this point it is worth noting that we have transformed a stochastic coupled map lattice, that 
involves a threshold condition and presents a clear separation of time scales, into a nonlinear stochastic 
partial differential equation, where the nonlinearity of the threshold is described by an infinite series 
of powers and the randomness enters via a Gaussian process, with zero mean to account for the 
dissipation at the boundaries. Along this transformation, and due to the approximations we have 
performed concerning the noise correlation, we have broken the time scale separation since the noise acts 
constantly in time. Nevertheless, we expect that such an equation explains the dynamical properties of 
the system within the fast time scale of the propagation of the avalanches. As we have mention before 
and it is discussed in |3l| , to deal with the slow time scale, where the avalanches are instantaneous, 
another noise correlation is more appropriate. 



III. DYNAMIC RENORMALIZATION GROUP PROCEDURE 

The model to be studied by the DRG is defined by the nonlinear partial differential equation ([l3]) and 
the Gaussian noise given by (^) and (^|). As a previous step we can check the relevance of the different 
coupling constants in this equation by naive dimensional analysis: a change of scale b = e > 1: 

P = e~ l f, t' = e- zl t, E' = e~ xl E, (16) 

is performed in Eqs. ( |l3| ) and (Q), being \ t ne roughness exponent, which is related to the hydro- 
dynamic exponent, and z the dynamical exponent. Then one obtains that the parameters transform 
as 

D^b z ~ 2 D, r ^b 2< - z -^- d T, X n ^b z+( -"- 1 '> x ~ 2 X n . (17) 

Under this scaling transformation, z and x are chosen to keep the linear model scale invariant, i.e., 
the parameters D and T have not to be modified. This choice gives z = 2, \ — (4 — d)/2, and 

\ n ^b^-V\ n . (18) 

Then one can see that when we apply iteratively the transformation (b — > oo) for d > 4 all the nonlinear 
terms vanish and they are irrelevant. However, all the coupling constants go to infinity for d < 4, and 
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hence all nonlinear terms become relevant; this implies that the upper critical dimension is d c = 4, 
and nontrivial values of the exponents are expected below it. 

The relevance of all the terms makes our problem much more complicated than for instance the 
Kardar-Parisi-Zhang model of interface growth, where only the first nonlinear term is relevant p2[ . 
The appropriate treatment of Eq. (|l^) would be to renormalize the infinite number of relevant coupling 
constants that are involved. Of course this is impossible to do in practice. In 0] an expansion in the 
number of coupling constants for the BTW model was performed with only odd terms, i.e., without 
symmetry breaking (K = 0). The critical exponents where obtained as a function of the highest 
coupling constant, up to Ag. Fortunately, the dynamical exponent was well behaved and could be 
extrapolated up to Aoq. However, keeping the symmetry of the step function, the nontrivial fixed point 
of the DRG is unreachable using the parameters given by Eqs. @ and @, even for the Zhang's 
model. We want to show that with the proposed alternative regularization of the step function, that 
breaks the symmetry of the BTW model and allows the existence of even coupling constants, the DRG 
fixed points are not changed but now they are accessible to the flow when the parameters take their 
real values. For this reason, and as a first attempt to justify our hypothesis as well as the conclusions 
of Ref. j24|, we will focus on Eq. ( [l3|) with only its first two nonlinear terms, i.e., A2 and A3, and see 
how they behave under a DRG transformation, 

^A^DV 2 E(r,t) + X 2 V 2 E 2 (f,t) + X 3 V 2 E 3 (r,t) + r 1 (f,t). (19) 
ot 

The DRG procedure consists of the removal of the fast modes (large wavenumber k) in the momentum 
space, followed by a rescaling of a factor e l in order to recover the original Brillouin zone [p5|-p7| . After 
this transformation, one obtains an equation that is equivalent to the original one but with different 
(effective or renormalized) coefficients. Successive iterations of this transformation give the flow of 
the coefficients in the parameter space. If this flow converges towards a fixed point, the system 
presents "scale invariance" in the hydrodynamic limit (large-distance and long-time behavior). Then, 
the fluctuations of the order parameter verify the scaling equation 

< [E{f ,t )- E{r + f,t + t)} 2 >^ r XF(t/r z ), (20) 

where the critical exponents \ an d z are those that ensure the existence of the fixed point. However, 
it is worth mentioning that with this procedure the scaling function F(x) remains unknown j33|. 
We now outline the DRG calculation. First of all we write Eq. ( |l9| ) in Fourier space: 

E = GoV- Goj^E *E-G E * E * E. (21) 

Here E(k,uj) and ry(fc,w) are defined as the Fourier transforms of E(f,t) and 7/(r, t), i.e., 

E(k,u) = J d d r dt e^*'^ E(f, t), (22) 

whereas 

Go(k,w)= . 1 (23) 
— iu> + Dk z 

is called the bare propagator. The symbol "*" represents the convolution product, defined as 

(E*E)(k,uj) = J d d qdfl E(q,Q)E(k-q,ui-fl). (24) 

Fig. 1 (a) shows the expression of Eq. (|2l]) in terms of Feynman diagrams. As the intensity of the 
noise T has also to be renormalized by the DRG transformation, we need to consider the equation for 
the correlation function of the transformed energy < E{k,ui)E{k' >, which, up to one-loop order, 
is: 

< EE' >= G Q G' a < W > + X f^ 2 2 ^'° <(E*E)(E* E)> >, (25) 



G 



where the prime denotes a dependence on k 1 ,u>' instead of the dependence on k, ui. The diagrammatic 
representation of this equation is shown in Fig. 1(b). Eqs. ( pl| ) and (|25|), that are the ones that we 
are going to renormalize, hold for < k < A, where A is the wavenumber cutoff due to the underlying 
discrete structure. The transformed noise rj(k, uj) turns out to be also a Gaussian process with zero 
mean, but with correlation 

< r](k, u>)r){k', J) >= 2{2Tt) d+2 T5 d {k + k')5(uj)S(cj'). (26) 




Fig. 1: A. Corral and A. Diaz-Guilera, "Symmetries and Fixed Point Stability of Stochastic Differential Equations Modeling 
Self-Organized Criticality" , submitted. 



FIG. 1. Diagrammatic expressions for Eqs. ( |2l| ) and (|2jJ), defined in the range < k < A. The double bar 
with the cross x at its end is the order parameter E, the single bar with the cross represents Gof?, whereas the 
single bar alone is Go- A vertex with n branches (n = 2 or 3 in the figure) represents a convolution product of 
n elements, including a prefactor — A n fc 2 /(27r)' n_1 '^ £i+1 - ) . The circles correspond to the average over the noise. 
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FIG. 2. Diagrams obtained after the first step of the DRG transformation. Now continuous lines correspond 
to the inner shell, whereas dashed lines correspond to the outer shell. The comparison with Fig. 1 allows to 
define new coefficients. Observe that the new averages only affect the outer shell. The notation has been 
simplified in respect to Fig. f , suppressing the symbol x at the end of the vertices and also the arrows. 

The first step of the DRG transformation consists in splitting the Fourier space in two shells: an 
inner shell, that contains the slow modes, i.e., < k < e~ l A, and an outer shell, containing the fast 
modes, e~ l A < k < A. Both modes are coupled through the convolution products in ( |2l| ) and (25) 



We consider the diagrams for the slow modes and perform a perturbative expansion of the fast ones 
up to the lowest order in the intensity of the noise (see Appendix for more details). Then we integrate 
out these modes by an average over the noise in the outer shell. After this transformation the resultant 
equations are shown diagrammatically in Fig. 2. It is clear that we can obtain new equations which are 
formally equivalent to the initial ones, ( ^l|) and (EjJ), defining the new coefficients as the original ones 
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plus the corresponding integrals over the outer shell. With the noise correlation (26) these integrals 
can be easily computed in the hydrodynamic limit (k — > 0, lo — > 0), as it is shown in the Appendix, 
and then the coefficients transform according to 



r. 



(27a) 



(27b) 



Ao — ► A2 



1 - 18 



£> 3 



12- 



L» 4 



(27c) 



1 - 18 



i d r\ 3 

D 3 



72 



D 4 



32 



D^X, 



(27d) 



where 



Id(l) = 



2S d l- e -'( d - 4 ) 
(2n) d d-4 



A 



d-4 



(28) 



and S d is the complete solid angle in d dimensions. However, the new equations are only defined in 
the inner shell < k < e~ l A. The second step allows to recover the original Brillouin zone < k < A 
by rescaling the equations using transformation (|l^) , which in Fourier space writes 



= e l k, 



E' = e -(x +z+d »E. 



(29) 



The combined effect of both transformations, in the limit / — > 0, constitutes an infinitesimal DRG 
transformation, which gives the flow equations of the parameters in parameter space. In these flow 
equations instead of A2 and A3 it is suitable to use the dimensionless coupling constants A2 and A3, 
given by 



Xi = 



7fTA| 
£ 4 



and 



A ^ 



(30) 



where = (dl d /dl) l=0 = [2S d /(2Tr) d ]A d - 4 . Then 



dT 

~dJ 



r [2z -2 X -d], 



(31a) 



dD 

~dT 



D\z-2-AXi 



3A 3 



(31b) 



dM 
dl 



= A, 



4- d 



20X\ 



24A, 



(31c) 



dAg 
dl 



= A, 



4-d- 



MXl 



Xf 

27A 3 - 32^ 
A3. 



(31d) 



We are interested in the invariance of the parameters under DRG transformations. This means that 
we have to look for the fixed points of the flow equations; if we write Eqs. ( |3l| ) as dcnj dl — gi(...otj ...), 
where ctj represents any coefficient, then the fixed points verify gi(...a*...) = 0, Vi. Considering D =/= 
and r ^ we obtain four algebraic equations with four unknowns, x, z, A2, and A3; their solutions 
will give us the fixed points of the transformation, A2 and A3 , as well as the values of the exponents 
z and x that guarantee that the DRG transformation leads to a scale-free behavior. Notice that the 
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particular values of T and D play no role in the existence and location of the fixed points. We can also 
find the stability of the fixed points under small perturbations using a linear stability analysis: the 
fixed point {a*} is stable (i.e., an attractor) if all the eigenvalues (or their real parts) of the matrix 
dgi/dctj evaluated at this fixed point are negative. 

The results are the following: for d > 4 one obtains six different fixed points, but the only stable 
one corresponds to 

X =~ z = 2, A* = A*=0, (32) 

where as usual e is defined as e = d c — d = 4 — d. This is the trivial or Gaussian fixed point, that 
gives a normal (or Brownian) diffusive behavior because of the vanishing of the coupling constants. 
The values of the exponents do not correspond with those of the Edwards- Wilkinson model, used in 
the study of surface growth, because the noise correlation is different p4j . For d < 4 this fixed point 
becomes unstable and the only stable one is 

18 9 2 3 27 

that was unstable for d > 4. In this case the diffusion is anomalous, to be more precise, the fact that 
z < 2 gives a superdiffusive behavior in the hydrodynamic limit. Note that the one-loop expansion 
in the intensity of the noise T gives a nontrivial fixed point which is expressed as a perturbation of 
the Gaussian one in a first order e-expansion. Observe also that the breaking of symmetry does not 
modify the value of the fixed point obtained without taking into account the even coupling constant 
A2 |Q. Moreover, the fact that the nontrivial fixed point is an attractor of the dynamics contrasts 
with equilibrium critical phenomena, where this point is stable only along one direction. In this fact 
lies the difference between fine tuning of parameters for equilibrium systems at the critical point and 
self-organization towards criticality for nonequilibrium processes. 

Now we know the attractors in the parameter space, but this is not enough in our case; since our 
stochastic equation (jl^) is derived directly from the discrete rules of the BTW and Zhang's models, we 
also need to know the basins of attraction of the stable fixed points and whether our initial conditions, 
that is, the initial values of the coefficients corresponding to our physical problem, are inside these 



basins. These values for the dimensionless coupling constants (30) can be calculated from Eqs. ( J14| ) 
and (Ua), and they are 



i /Wr°/ (2 W - li^f^ixi 

- 4 a 2 EZ fW(K)*' 3 6 a 2 £ c 2 /«(f) 3 ' 



result that also holds for K = 0, where we obtain AJJ = even for the Zhang's model. The superscript 
"o" indicates the initial value of the coefficient, that is, its value before any renormalization. As we 
have no restriction for r° (except that it has to be small), and K can take any arbitrary real value, 
this implies that the initial dimensionless coupling constants will be defined in the following region: 

^ < ^(A^) 2 , (35) 

having used for f(x) the explicit form given by Eq. (|^). 

Clearly, the stable fixed point for d < 4 ( |33| ) is outside the region of initial conditions defined by 
Eq. (^). It will be of the maximum interest however to know whether these conditions will drive 
the system towards the nontrivial fixed point or not. We first consider K = 0, which implies ASJ = 0, 
corresponding to the case studied in Ref. ^4|. For d < 4 one gets a different behavior depending on 
the sign of A3. Fig. 3 shows that when A3 is positive it flows towards the stable fixed point A3 = e/27 
giving a dynamical exponent z = 2 — e/9. A negative A3, which is our case of physical interest, flows 



away. An exact solution of Eq. (31d) with A2 = gives that A3 would reach —00 in a finite I and then 
would reappear as A3 = 00, being then under the attraction of the nontrivial fixed point. However, our 
one-loop calculation forces the flow of the coupling constant along the parameter space to stay of order 
e, and one cannot sustain the validity of the preceding description. Then it is not possible to predict 
the renormalization of A3. It will be either renormalized to A3 or other fixed points will appear along 
the flow (corresponding to strong coupling and not given by the one- loop e-expansion). Therefore, the 
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fixed point given by Eq. (33), although it is an attractor, is unreachable from our initial conditions 
(A3 < 0). For that reason the conclusions of Ref. [Q were incomplete. On the other hand, above the 
upper critical dimension the system evolves towards the trivial fixed point A3 = giving a diffusive 
behavior with z — 2 provided that A3 is not too much negative (see Fig. 3). This behavior of the fixed 
point A3 as a function of e corresponds to a transcritical bifurcation. 



u s 
-o »• o- 



u s 
-0 »■ o- 



Fig. 3: A. Corral and A. Diaz-Guilera, "Symmetries and Fixed Point Stability of Stochastic Differential Equations Modeling 
Self-Organized Criticality" , submitted. 



FIG. 3. Flow in A3 space when only this nonlinear term is taken into account. The squares correspond to 
the stable (S) and unstable (U) fixed points and the arrows show the flow under DRG transformations. 

Now, by introducing the alternative regularization (K ^ 0), we will see the effect of the symmetry 
breaking. First of all, we insist that the stable fixed points are the same as for K = 0, due to the fact 
that A2 renormalizes to zero. Moreover, as can be seen in Fig. 4, where we have plotted the flow lines 
of Eq. fey) obtained by numerical integration, the basin of attraction of the nontrivial fixed point is 
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delimited by the parabola 



A - 4 A 2 



(36) 



which is also a particular solution of the flow equations, no matter the value of d. This parabola is 
inside the region defined by Eq. (|3^) , and this fact implies that the new regularization makes possible 
to reach the attractor for d < 4 starting in the region of physical meaning. Using Eqs. |34|) and ( |36| ) 
together with (|9|) one gets that the condition to converge towards the nontrivial fixed point is K 2 >\- 
Then, the parameter that breaks the symmetry in the regularization of the step function, which in 
principle was arbitrary, determines the behavior of the system in the hydrodynamic limit. 



0.15 



0.05 




-0.05 



A 2 

FIG. 4. Flow in (A2, A3) space for d — 3 when both nonlinear terms are taken into account. In gener al fo r 



the same. Dots correspond to the numerical integration of Eqs. (31c) 
!q), that clearly delimits the basin of attraction of the nontrivial fixed 



any d < 4 the results are qualitatively 
and (31d), and the thin line is Eq. ( 
point, as it is seen in theplot. Below the continuous thick line the values of the parameters correspond to our 
physical situation, Eq. (p5|). Squares correspond to the fixed points. Observe that for A2 = we obtain the 
same results as in Fig. 3. 



For d > 4 the flow is more complex, because of the six fixed points, but the result is that convergence 
towards the Gaussian one also happens for our initial conditions, as Fig. 5 shows. The linear stability 
analysis of the fixed points gives the same results as the numerical integration shown in the figure. 
However this linear analysis fails for d = 4, where all the fixed points collapse towards the Gaussian 
one. It is by means of the numerical integration that we verify that it is an attractor for the region 
above the parabola given by Eq. (|35|), but for the region below it is a repeller. This strange behavior 
appears because in d — 4 we are at the bifurcation point. 
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0.1 0.2 0.3 0.4 0.5 0.6 0.7 



A-2 

FIG. 5. Same as Fig. 4 but for d = 5. The results hold for d > 4. Only four of the six fixed points are 
shown because of the symmetry of the flow lines. In this case, the curve A3 = — §A|, represented by another 
thin line, is the repulsive branch of the saddle point. 

In Ref. |24|] it was shown for the BTW model and K = that A2n = 0, whereas for the Zhang's 
model, although the even coupling constants do not vanish, it was argued that their flow equations 
became decoupled from the odd ones in the limit (3 — > 00. This fact enabled to establish the same 
universality class for both models, and to deal with only odd terms in Eq. (|ll|). Then, an expansion in 
the number of coupling constants was performed, whose extrapolation compares well with the results of 
the simulations . Note that in the simulations one computes the dynamical exponent relating the 
characteristic length and lifetime of the avalanches, whereas within the DRG framework one computes 
the dynamical exponent from the fluctuations of the order parameter p5| . The agreement between 
these calculations confirms the basic scaling hypothesis that in both cases length and time are related 
by means of the same exponent. However, the problem of this calculation was that the nontrivial fixed 
point was unreachable for the original equation. 

In our approach, due to the symmetry breaking, we have to consider also the effect of even coupling 
constants. In the present work we have dealt with a restricted problem with only the lower-order even 
and odd coupling constants, A2 and A3, showing that A2 renormalizes to zero, supporting the calculation 
of Ref. [^4| . Then the stable fixed points are not modified by the presence of an even coupling constant 
in the model, but due to the symmetry breaking that we have introduce, the nontrivial one is an 
attractor in the parameter space when the parameters corresponding to the real model are taken into 
account. This behavior should be the same for any even coupling constant; actually, preliminary 
calculations including A4 and A5 in Eq. ( JToj ) make us to suspect that all even coupling constants 
renormalize to zero. This fact means that in the hydrodynamic limit the solution of both models has 
to be symmetric under parity transformations of the order parameter; then, for the BTW model the 
DRG restores the broken symmetry, whereas for the Zhang's model we conclude that its asymmetric 
nature is irrelevant in the behavior at large distances and long times. Therefore this validates the 
extrapolation performed in pi| since now we have showed that the symmetry breaking makes the stable 
fixed points reachable, when starting in the region of physical interest in the space of parameters. Let 
us finally mention that in a recent work, Ghaffari and Jensen J3(| perform a different extrapolation 
of the same results which show a better agreement with large-scale simulations and with real-space 
renormalization calculations for the dynamical exponent [ ^3| . It is noticeable than the same technique 
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has been applied to the study of the effect of dissipation in a uniformly driven BTW model 



IV. CONCLUSIONS 

We have studied analytically two models which show self-organized criticality. The difference between 
them is that the second one (BTW) is symmetric under a parity transformation, whereas the first 
(Zhang's) model is not. From the microscopic rules one writes a effective long wave-length equation 
involving the threshold condition, which enters into the equation through a step function, making 
the equation unapproachable under this form. We have introduced a new regularization of the step 
function that breaks the symmetry of the BTW model. After a power series expansion, the equation 
is suitable for the application of the dynamic renormalization group, although it contains an infinite 
number of relevant coupling constants. In consequence one has to truncate at some point the expansion 
in the coupling constants. The results only have sense if it is possible to extrapolate the values of 
the exponents up to an infinite number of coupling constants. We obtain the fixed points of the 
transformation in parameter space and study carefully their stability and basins of attraction. Then 
we find that with this regularization it is possible to reach the nontrivial fixed point for d < 4, that was 
unreachable in a previous work, where symmetry was not broken. This means that in the hydrodynamic 
limit the models display "scale invariance". Moreover, in this limit we obtain a symmetric behavior 
under parity transformations for both models, and therefore the recovery of the broken symmetry for 
the BTW model and the irrelevance of this symmetry for Zhang's one. Although we have dealt with 
a simplified version of the problem, we expect this behavior to be the same for the complete problem, 
in the sense that all even coupling constants renormalize to zero, validating the calculation of Ref. 
[ p4[ . The application of this technique should also be useful for other kinds of problems in which one 
deals with thresholds or with an infinite number of nonlinear terms, for instance interface dynamics. 
Moreover, the performed DRG calculation is interesting because it provides an example showing how 
much important is to know not only the stable fixed points of a DRG transformation but also their 
basins of attraction. It is remarkable the fact that a simple symmetry breaking can solve the problem 
of the inaccessibility of the attractors in parameter space. 

The fact that the parameter that breaks the symmetry determines the behavior in the hydrodynamic 
limit is difficult to understand and we believe that it is an artificiality introduced in the calculation by 
the truncation in the coupling-constants expansion. We expect that higher orders in this expansion 
will give a behavior independent on the K value. 
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APPENDIX: 

Here we present further details about the derivation of Eqs. (^7|), which give the transformation of 
the parameters after the first step of the DRG. Our starting points are Eqs. ( ^l| ) and (|25|), i.e., the 
equations for E(k,uj) and < E(k,uj)E(k' , u/) >. As we have already mentioned, these equations are 
only defined for < k < A. The DRG procedure consists of splitting the momentum space into two 
shells, an inner one, with < k < Ae~ l , and an outer one, with Ae~ l < k < A. Then the magnitudes 
that depend on fc, like the energy E, split in the following way 

E{k,oj) = E < (k,u)+E > (k,oj) = E{k,uj)Q{Ae- 1 - k) + E{k,uj)Q{k - Ae~ l ) (Al) 

where Q(x) is again the Heaviside step function. This equality defines E < (k,u>) as the corresponding 
part of the energy in the inner shell, whereas E > (k,uj) is the same but defined in the outer shell. This 
separation also holds for the bare propagator Go and the noise rj. 
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The DRG procedure eliminates the modes of the outer shell, within the same philosophy than 
the Kadanoff transformation does in real space. Then one is only interested in E < (k,uj) and 
< E < (k ) u>)E < (k',oj') >, whose equations turn out to be equivalent to j2l| ) and but with 

additional terms due to the coupling between the two shells, via the convolution products. The 
fact that E > (k,u>) appears in the inner-shell equations allows a perturbative expansion in the form 
E > (k,uj) — Gq (k, u))rj > (k, us) + ... (using the equivalent of Eq. (|2l| ) but in the outer shell). Then 
the noise in the outer shell enters into the equation for E < (k,u). A similar perturbative expansion is 
done for < E < (k,u})E < (k' ,u>') >. By averaging over 7y > (A;, u>), the contribution of the fast modes is 
eliminated from the inner shell. This is done up to one-loop order in the perturbative expansion, that 
is, the lowest order in the intensity of the noise T, which implies that it has to be small enough. This 
tedious calculation becomes more appealing using the diagrams of Fig. 1 , instead of the corresponding 
equations. After this process, the relevant diagrams that survive the averaging are shown in Fig. 2. 
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Fig. 6: A. Corral and A. Diaz-Guilera, "Symmetries and Fixed Point Stability of Stochastic Differential Equations Modeling 
Self- Organized Critically", submitted. 



FIG. 6. Diagram computed in the Appendix as an example. The brackets stand for an average over the 
outer shell. 



As an example let us consider one of them, shown in Fig. 6 and denoted by V(k,cu): 

V(k,u) =< ^^G<(k,cj) J d d q dQG>(k-q,LJ-Q)G>(q,n)rj>(q,n) x 
-A 2 (fc-<7) 2 f , d , ,_, 



(27T) 



d+1 



d d q' d£VG$ (k - q- q' ,u - fi - Q.')E < {q' x 



" A2 (2tt)^i J ddq " dn " G ° (k-q-q'-q 7 '^-n-n'- fi") x 
X77>(fc - q- q> - q>>,LU - Q - Q' - Q")E < (q 7 ', Ct") >>, 



(A2) 
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where the symbol < ... >> stands for an average over the outer shell. Using the noise correlation given 
by Eq. ( p6| ) we can integrate over f2, Ct", and q", and then we have 

2A 3 r f - 

V{k,uj) = - {2ix y d+1 k 2 G<{k,Lo) J d d q' dQ'E<(q',n')E<(k-q',LJ-n')x 

x / d d q G>(k - q,u,)G>(q,0)$ - q) 2 G>(k - q- q',oj - fi')(fc - q- q') 2 G>(q,0). (A3) 



As the bare propagator is a known function, given by Eq. (p3|), we are also able to perform the integral 
over q, that is 



d d q 



d d q G>(k - q,co)G> Z (q,0)(k - q) 2 G>(k - q- q',u- fi')(fc - q- q>) 



(A4) 



This integral is a function of k, ui, q', and Q' . However, we are going to evaluate it in the hydrodynamic 
limit, by taking k, q' — * and u>, il 1 — > 0. Then 



d d q 



d d q q 4 G> 2 (-q,0)G> 2 (q,0) 



,d-4 



— [ Q d - 5 da=—- (l- e - l{d -^ 

D*J A e-< D*d-4\ 



(A5) 



It is easy to check that this result is also valid for d = 4. We have used the explicit form of the bare 
propagator (|2^ ) and also that d d q = S c iq d ~ 1 dq, with Sd the complete solid angle in d dimensions, that 
is, the area of a unit (d)-sphere. Then, by making use of Eq. we obtain 



d d q 



(2n) d I d (l) 
2D 4 '■ 



and substituting into Eq. (AS), 
V(k -> 0,lj -> 0) 



.—Ar^G^(k,u){E<*E<){k,u)I d (l)^ r 



(A6) 



(A7) 



It is clear from Fig. 2 that after the first step of the DRG we have the same as at the beginning 
(Fig. 1), but defined only in the inner shell, plus a lot of diagrams of the same type as the one in Fig. 
6. These diagrams, which contain integrals over the outer shell, renormalize the other diagrams that 
are only defined in the inner shell. For instance, if we consider the following diagram 




(A8) 



and compare it with with Eq. (A7) we observe only an additional term I^tyX^r / D 4 that comes from 
the outer shell integration. So, the diagram shown in Fig. 6 contributes to the renormalization of 
([AS]), that is, renormalizes the coupling constant A2. As Fig 2(a) states, the diagram in Fig. 6 appears 
eight times in the perturbative expansion, and the new A2 , after the first step of the transformation 
will be modified by 



A. - A.( 1 + 8J„(I)M£ 



(A9) 



In the same way one can perform the outer-shell integrals of the rest of diagrams in Fig. 2(a). A 
general result for its contribution to the renormalization of any coupling constant A„ or to the diffusion 
coefficient D (which will be referred here also as — Ai) is given by 
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(-ir-^G) 



r IT 



D v+1 



A, 



(AlO) 



where v is the number of vertices each diagram has, three for our example (since the dashed line in 
Fig. 6 forms a triangle), b(m) is the number of branches of the m-th vertex, two for each one in the 
example, and B is the number of branches of the diagram that is renormalized (two in the example), 
and fulfills B = J2m=i b{m) — v — 1. Note that th e magnit ude in Eq. ( AlO ) is dimensionless. Using 
this equation and Fig. 2(a) the derivation of Eqs. ( 27b - 27d ) is then straightforward. 

For the renormalization of the intensity of the noise T we have only one diagram, the dashed one in 
Fig. 2(b). It is immediate to see that the integral over the outer shell, no matter is value, is multiplied 
by a factor k 2 k' 2 . Then 



r l 



A\?k> 2 % 



(AH) 



where A is simply a numeric factor, and hence, in the hydrodynamic limit and up to one-loop order, 
the intensity of the noise is not renormalized, after the first step of the DRG, as Eq. (27a) states. 
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